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The effects of uniform horizontal shear on a stably stratified layer of gas is studied. The system is initially destabilized by 
a magnetically buoyant flux tube pointing in the cross-stream direction. The shear amplifies the initial field to Lundquist 
numbers of about 200-400, but then its value drops to about 100-300, depending on the value of the sub-adiabatic gradient. 
The larger values correspond to cases where the stratification is strongly stable and nearly isothermal. At the end of the 
runs the magnetic field is nearly axisymmetric, i.e. uniform in the streamwise direction. In view of Cowling's theorem the 
sustainment of the field remains a puzzle and may be due to subtle numerical effects that have not yet been identified in 
detail. In the final state the strength of the magnetic field decreases with height in such a way that the field is expected to 
be unstable. Low amplitude oscillations are seen in the vertical velocity even at late times, suggesting that they might be 
persistent. 
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1 Introduction 

Dynamos convert kinetic energy into magnetic energy. A 
typical example is thermally driven convection. If the dif- 
ference between heating and cooling across the domain is 
strong enough, thermal energy can be converted into kinetic 
energy by the Rayleigh-Benard instability, and part of this 
kinetic energy can then be converted further into magnetic 
energy by the dynamo instability. However, if there is shear, 
then this can supply a major part of the kinetic energy. In 
view of Cowling's (1933) anti-dynamo theorem, it is clear 
that the magnetic field must be fully three-dimensional. This 
is easily achieved if the motions are three-dimensional as 
well. A linear shear flow is just one-dimensional, although 
it is in principle possible that this flow becomes nonlinearly 
unstable and develops fully three-dimensional turbulence. 
Unfortunately, this requires rather large Reynolds numbers 
and is not easy to achieve. However, in a stratified system it 
is also possible to produce three-dimensional motions by the 
magnetic buoyancy instability. This was an important agent 
responsible for driving a dynamo in the model of Cline et 
al. (2003). Yet another possibility is the magneto-rotational 
instability, which can also produce three-dimensional mo- 
tions to drive a dynamo (Brandenburg et al. 1995; Hawley 
et al. 1996; Stone et al. 1996), but this requires that there is 
also rotation. 

In this paper we investigate a system similar to that of 
Cline et al. (2003), but with a linear shear profile instead of 
a sinusoidal one, or rather a modified sinusoidal one where 
one flank is steeper than the other Another difference is that 
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in our model the shear extends through all layers and is not 
localized near the bottom of the domain, as in the model of 
Cline et al. (2003). Their setup was motivated by the pres- 
ence of a strong shear layer in the solar tachocline (Hughes 
et al. 2007). However, here we are interested in more gen- 
eral aspects rather than particular applications. It is possible 
that some important dynamics would be lost by ignoring 
the additional vertical dependence of shear, but this is not 
well understood at present. Regardless of whether or not 
self-sustained dynamo action exists, there are a number of 
issues that deserve to be addressed in such a setup. Firstly, it 
is useful to determine the energy fluxes between kinetic and 
magnetic energies and how they are coupled to the shear, 
for example. Secondly, if magnetic buoyancy plays a role in 
producing turbulence, it should be possible to quantify this 
by measuring suitable correlations. In particular, it is not 
clear whether a stratification that is close to marginally sta- 
ble is advantageous compared to one that is strongly stable. 
Finally, in order to assess the possibility of large-scale dy- 
namo action, one needs to determine the turbulent transport 
coefficients. 



2 The model 

2.1 Governing equations 

We consider a Cartesian domain with externally imposed 
linear shear and vertical gravity leading to density stratifica- 
tion in the z direction. The full set of hydromagnetic equa- 
tions for the magnetic vector potential A, the density p, the 
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velocity U, and the specific internal energy e, can then be 
written in the form 

VA 

'Vt 

Vlnp 



Vt 



= -SAyX - {VUfA - fiovJ, 
= -V-U, 



(1) 
(2) 



VU 

p—=-SpU^y-Vp + pg + JxB + V-2upS, (3) 

P^ = -pV • [/ + V • KVT + 2pvS'' + txovJ^ , (4) 

where V/Vt = d/dt + {U + Uq) ■ V is the advective 
derivative, Uq = (0, Sx, 0) is the imposed large-scale shear 
flow, Sij — ^{Uij + Uj^i) — ^(5ijV • C is the traceless 
rate of strain tensor, B = V x A is the magnetic field, 
J = V X B//io is the current density, is the magnetic 
permeability, t] and v are respectively the magnetic diffusiv- 
ity and kinematic viscosity, K is the heat conductivity, and 
g = —gz is the gravitational acceleration. The fluid obeys 
an ideal gas law p = pe(7 — 1), where p is the pressure, and 
7 = Cp/cv = 5/3 is the ratio of specific heats at constant 
pressure and volume, respectively. The internal energy per 
unit mass is related to the temperature via e = c^T. 



2.2 Initial and boundary conditions 

Our initial stratification is a polytrope where p is propor- 
tional to p^, and r is related to the polytropic index m via 
r = 1 + 1/m. The superadiabatic gradient is usually defined 
as the normalized entropy gradient with respect to logarith- 
mic pressure, V — Vad = d(s/cp)/dlnp, i.e. 



V- V, 



ad 



7 



(5) 



A stable stratification corresponds to V — Vad < 0, i.e. 
r < 7. Following Cline et al. (2003) we adopt m = 1.6, 
which gives V - Vad » -0.015. (We note that the V - Vad 
quoted by Chne et al. (2003) is scaled by a factor 2(m -|- 1), 
giving -0.08.) 

In a polytrope, the temperature, and hence the square 
of the sound speed are proportional to the negative gravity 
potential. Using = (7— l)/i we find 

.2 7 



m+1 



(6) 



where $ = —{z^ — z)g is the gravitational potential, and 
the Zoo > 2; is the top of the atmosphere where temperature, 
density and pressure would vanish. In order that p = po and 
Cg = CgQ at a certain reference height z = —Ho,-we choose 



= -Ha 



19 



(7) 



We consider a computational domain with horizontal extent 
-2 < x/Ho < 2,-2 < y/Ho < 2, and vertical extent 
—4.12 < z/Hq < —1, where z/Hq = —1 corresponds to 
the top of the layer. 



Our initial magnetic field is given by a flux tube pointing 
in the x direction, but with a perturbation in the x direction. 
Thus, our vector potential is given by 



. „ 1 + e cos kxX , 
A = Bq ^ . n , „o X X or 



(8) 



1 + dr-^/E? 

where 5r = [x, y, z — zq) is the distance from the core of 
the tube in the yz plane. We choose Bq = Q.I^qHqpqpq, 
e = 0.3, R = 0.2iJo, and zq = -3Hq. 

In all cases we use stress-free boundary conditions for 
the velocity. 



0, 



together with a vertical field condition, i.e. 



B^, = B,, 



0. 



(9) 



(10) 



For the specific entropy we use either an extrapolating 
boundary condition, which allows the values of thermo- 
dynamic variables on the two boundaries to change freely 
(Run A), or we fix the energy influx at the bottom and the 
temperature at the top (Run B). The latter condition is also 
used by Cline et al. (2003). 

2.3 Units and control parameters 

Non-dimensional quantities are obtained by setting 

d = 9 = Po = Cp= IJ,Q = 1, (11) 

where po is the initial density at Ztop- The units of length, 
time, velocity, density, entropy, and magnetic field are 



[x] = Hq , [t] = ^ fH^g, \U \ = \p\ = PO , 

[s] = Cp , [B] = ^JgHopQiiQ. 



We define the fluid and magnetic Prandtl numbers as 

V „ V 



Pr — P 



77 



(12) 



(13) 



where xo = K/[p^Cp) is the thermal diffusivity, and po is 
the density at ^; = —Hq. We define the magnetic Reynolds 
number and the shear parameter via 



Rm = 



Sh = 



(14) 



where fcf = 2it/Hq is assumed as a reasonable estimate for 
the wavenumber of the energy-carrying eddies. 

The simulations were performed with the PENCIL 
CODE^ which uses sixth-order exphcit finite differences in 
space and third order accurate time stepping method. The 
dependent variables are A, In p, U , and the specific entropy 
s, which is related to the speed of sound Cg via 



[(7 - 1) ln(p/po) + 7s/( 



•pj) 



(15) 



where Cgo and po are normaUzation constants. 



' http://www.nordita.org/software/pencil-code/ 
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Fig. 1 Final stratification for Runs A (solid lines) and B (dashed 
lines). Note that Run A is nearly isothermal (constant sound speed 
Cs) while Run B is nearly isentropic (V — Vad is close to zero). 

3 General features of the simulation results 

In the following we discuss Runs A and B that differ only in 
the boundary conditions adopted for the entropy. For Run A 
we use extrapolating boundary conditions for the specific 
entropy, while for Run B we fix the energy influx at the 
bottom and the temperature at the top. In all cases we use 
the same polytropic initial condition. 

We have evolved both simulations for several thermal 
and magnetic diffusion times. The mean stratification set- 
tles then to a new stratification (Fig. 1). Runs A and B differ 
mainly in the mean temperature gradient which becomes 
nearly constant in Run A and stays finite for Run B. For 
Run B we find that V — Vad ~ —0.01; see Fig. 1, so it is 
stably stratified, but close to marginal. By contrast, in Run B 
the stratification is strongly stable with V — Vad ~ —0.38. 

When the simulation is started the initial flux tube be- 
gins to rise owing to magnetic buoyancy, and it is also be- 
ing sheared out in the streamwise direction. Both effects can 
clearly be seen in a tz diagram of and By; see Fig. 2. 
Unless noted otherwise, an overbar denotes averaging over 
the X and y directions. The initial phase of the buoyant rise 
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Fig. 2 (online colour at: www.an-journal.org) Visualization of 
Bx and By as functions of t and z for Run A. Note the initial 
ascent of Bx as well as the subsequent amplification of By. For 
comparison, the curve z = ^{5p/ p)gt^ — zo has been overplotted 
with an estimated value of 5p/ p — 0.2. 

follows approximately a parabolic trajectory, as is indicated 
by the superimposed dashed line in the first panel. 

The effects of magnetic buoyancy can also be estab- 
lished at later times through a systematic correlation be- 
tween strong fields and upward motion. We define a nor- 
malized buoyancy parameter, 

Bu = (ipU,yB^)/{poU,^sBf^,), (16) 

where the dash in [pUz)' denotes the departure of the mass 
flux from the horizontal mean. It turns out that Bu is al- 
ways positive. This means that strong fields are systemat- 
ically correlated with upward motions and vice versa. For 
Run A this value is « 0.035 and for Run B it is w 0.07. This 
supports the idea that magnetic buoyancy is indeed active, 
but it is not clear that it still plays an important role at late 
times compared to early times. The peak values of Bu are 
around 3.5 at t sa b{g/Ho)^/'^ both for Runs A and B. 

In Figs. 3 and 4 we show the values of Uz and Bz on the 
periphery of the computational domain for Runs A and B at 
different times. Toward the end of the simulation the state 
appears nearly perfectly steady, but there are actually per- 
sistent low amplitude oscillations that remain excited at all 
times. In Fig. 5 we show a zt diagram of the mean vertical 
velocity, U z- The frequency of these oscillations is about 
0.9 (g/Hoy^^, This value lies between that of the Brunt- 
Vaisala frequency, whose local value varies between 0.4 and 
0.5 from bottom to top of the domain, and the acoustic fre- 
quency of about 2 {g/HoY^^. Alfven oscillations are also 
possible; their frequency is about 0.2 for Run A in the final 
state. 

In Fig. 6 we plot the evolution of C/rms and i?rms 
in diffusive time units, expressed in terms of the mag- 
netic Reynolds number R^^ = C^rms/?7^i and Lundquist 
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Fig. 3 (online colour at: www.an-joumal.org) Visualization of Uz (upper row) and (lower row) for Run A at different times. The 
coordinate directions are indicated in the upper left panel. Note the symmetry of and the antisymmetry of with respect to a:: = 0. 




Fig. 4 (online colour at: www.an-joumal.org) Same as Fig. 3, but for Run B. Here both Uz and B^ have no symmetry property with 
respect to a; = 0. 
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Fig. 5 (online colour at: www.an-joumal.org) Horizontally av- 
eraged mean vertical velocity showing persistent low amplitude 
oscillations for Run A at late times after t = 8200 {Ho/gY^'^. 
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Fig. 6 Comparison of the evolution of rms velocity and mag- 
netic field, normalized in resistive units and expressed in terms of 
Rm and Lu, for Runs A and B. The lower abscissa give time in 
dynamical units and the upper one in resistive units. 



number Lu = WA.rms/??^!, respectively. Here, "yA.rms = 
-Sims/ V/iopO is the Alfven speed. Note that toward the end 
of the simulation we have i?,n ~ 10 and Lu « 300. There 
is no evident tendency for decay. 



4 Flow of energy 

In order to understand where the energy comes from we 
consider here the energetics of the system. In a stratified 
hydromagnetic system there are four forms of energy: po- 
tential, kinetic, thermal, and magnetic energy, defined re- 
spectively by 



Ep= J p^dV, Ek = J ^pU\lV, 




-1 




Fig. 7 (online colour at: www.an-journal.org) Zonally averaged 
velocity (upper panel) and magnetic field vectors (lower panel) in 
the xz plane superimposed on the corresponding y components of 
these fields for Run A in the final and nearly steady state. Note the 
presence of a mean field with Bx > and By < 0, so B^By < 0, 
as expected for negative shear. 



Et = J pedV, Em^ j 2^B'dF. 



(18) 



In the presence of stress-free boundary conditions, the evo- 
lution of these energies is governed by the following four 
ordinary differential equations: 

Ep + W^b = 0, (19) 
Ek - M^b - VKc - Wl - e| + EK - 0, (20) 
+ VKc - eK - CM - ibot + itop = 0, (21) 

+ W^L - + EM - 0, (22) 
where dots denote time derivatives, 

W^b = JpU-gdV, Wc = JpV-UdV, (23) 

describe the work done by buoyancy and compression, 

Wl= Ju-{JxB) dV, (24) 

is the work done by the Lorentz force, 

4 = - / pU{SU)dV, 4 = y" Po'B{SB)dV, (25) 

couple shear to kinetic and magnetic energy reservoirs. 



(17) eK= 2piyS^dV, <J-^J^dV, 



(26) 
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Fig. 8 (online colour at: www.an-journal.org) Same as Fig. 7, 
but for Run B. 



are viscous and Joule dissipation, and 

itop = / i^R • dS, Lbot = / i^R • dS, (27) 

..'top J hot 

are the luminosities at top and bottom of the domain. The 
evolution of total energy, Etot = Ep + Ek + Et + Em, is 
therefore governed by 

s I 



L 



hot 



Lt 



top, (28) 

i.e. energy is supplied by shear coupling to kinetic and mag- 
netic energies, as well as by heating from below, and energy 
is removed by radiation at the top. 

Energy is exchanged and transformed into other types 
of energy, as is shown in Figs. 10 and 11 for Runs A and 
B, respectively. It should be noted that, even if the system 
is not perfectly in a steady state, we cannot expect perfect 
energy balance, because the code solves the equations in 
non-conservative form. However, the discrepancies should 
become smaller at larger numerical resolution. We notice 
that in all cases the energy flux at the top of the domain 
is less than at the bottom, and that these values are much 
larger than any of the other energy fluxes in the system. This 
is related to the discretization error of the scheme and does 
not indicate that energy is taken from the thermal energy 
flux. For Run B there is another obvious problem in that 
the magnetic energy dissipation is much larger than what 
can be accounted for by the influx of energy by shear. This 
is related to the fact that in this run there is a fairly strong 
magnetic field whose maintenance may be explained purely 
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Fig. 9 Plot of the energy input from shear via the terms ek ™d 
e^- Note the larger e? 
compared with Run B. 



e^- Note the larger excess of e^^ compared with for Run A 
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Fig. 10 Flow diagram showing the conversion between different 
energy forms for Run A. All energy fluxes are normalized with 
respect to the averaged magnetic energy content times the modulus 
of the shear rate and multiplied by a factor of 1000 to bring the 
values closer to unity. 



as the result of numerical errors, even though no obvious 
problems can be seen from images such as Figs. 3 and 4. 
This could be related to the fact that there are low ampli- 
tude azimuthal variations with a relative amplitude of about 
10^^ that cannot be seen unless the mean field is subtracted 
out. In any case, we should keep in mind that the reason for 
the maintenance of magnetic energy could be related to a 
very weak "non-axisymmetry" of the field that is related to 
numerical noise. Regardless of these subtle shortcomings, 
there are many aspects that do make sense. In particular in 
Run A the balance between inflow and outflow of magnetic 
energy is quite reasonable, and so is the balance of kinetic 
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Fig. 11 Same as Fig. 10, but for Run B. 

energy in Run B. This suggests that the flows in Run A are 
mainly driven by work done by the magnetic field against 
the shear, while in Run B the driving comes mostly from 
work done by the Reynolds stress against the shear. 

5 Profiles of mean field and velocity 

It is interesting to confront the mean profile of the magnetic 
field with stability criteria for magnetic buoyancy instabili- 
ties. Following Newcomb (1961), a necessary condition for 
instability is 

dln\B\/dz> jHpN^/v\ (29) 

(see also Hughes & Proctor 1988). This condition is actually 
met in the upper quarter of the simulation domain where 
there is a rapid decrease of |B| with height; see the top panel 
of Fig. 12 for Run A. 

The mean "toroidal" velocity, Uy shows a systematic 
variation in the z direction that is reminiscent of that found 
for isotropically forced shear-flow turbulence (Kapyla et al. 
2009), which was then interpreted in terms of a vorticity dy- 
namo (Elperin et al. 2003). In the present case, however, the 
magnetic field is already so strong that such an effect would 
be suppressed. Therefore, the flow is here more like a direct 
response to the magnetic field. The other two components 
of the mean flow are negligible by comparison. However, 
we recall that the z component of the mean velocity shows 
persistent low amplitude oscillations (Fig. 5) that are not 
present in the other two components. 

6 Mean-field transport coefficients 

In order to characterize the flow properties further we now 
consider the ability of the flow to mix and to produce large- 
scale magnetic fields. We do this by using the test-field 
method of Schrinner et al. (2005, 2007) with sinusoidal 
test fields, as explained in detail in Brandenburg (2005) and 




Fig. 12 By (solid line) together with 100 x Bx (upper panel) 
and Uy (solid line) together with 10 x Ux and 10 x Uz (lower 
panel), for Run A, time averaged over the last 1000 time units. 
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Fig. 13 (online colour at: www.an-joumal.org) Diagonal com- 
ponent of the a tensor for Run A at early times (upper panel is 
for axx and lower panel for ctyy). Note that ayy is negative in the 
lower part of the initial rising flux tube and positive above, while 
axx is negative throughout. 



Brandenburg et al. (2008a,b). In Fig. 13 we plot axx and 



^yy 



versus time and height for Run A at early times. It turns 
out that axx is mainly negative at the location of the initial 
rising flux tube while ayy is negative in the lower part of the 
tube and positive above. In order to interpret this result we 
compare now with expectations from mean-field theory. 

The combined presence of shear and magnetic buoyancy 
is particularly interesting, because one may expect there to 
be an additional contribution that is related to the magnetic 
stress, hxhy (Brandenburg 1998). Such behavior was alsore- 
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Fig. 14 Vertical profiles of —ThPuI^/g (solid line), 
bxhz,y/g{lJ,QpoY'^ (dashed line), and ebuoy x bxby/c^{fioPoy^^ 
with ebuoy =0.1 (dotted Une) for Run A at t = 10(i^o /gV^^- 

produced using the second order correlation approximation 
(Riidiger & Pipin 2000). In the spirit of the t approximation 
(e.g. Blackman & Field 2003) this term can be derived by 

calculating d£ / dt, which has contributions from u x b and 
ii X b. The first term, which gives rise to the usual kinetic 

a effect, is 

{U X b)y = ByU^-yU^ + (30) 

while the second term, which gives rise to the magnetic a 
effect as well as a new term proportional to b^by. To derive 
the second term we consider here the momentum equation 
in the form 

u = ... + B -b/piio + ^g, (31) 

where dots refer to additional terms that are less relevant for 
the present discussion. Replacing Sp/p by B^/2ppo and 
linearizing about By, i.e. 6p/p « —Byby/pjjbQ, yields 

(it X b)y = By (fo^,y6a;/p/io + S^'y^x/^op) + ••• • (32) 
Therefore the yy component of a should be 

^yy ^ -TUa;Uz,y-\ [bxbz.y + ebuoybxby/Ho)+...,(33) 

Pop 

where ebuoy is an empirical non-dimensional parameter 
quantifying the relative importance of magnetic buoyancy 
effects. In Fig. 14 we plot the contributions from all three 
terms for Run A at i = 10(ffo/5)^^^- It turns out that the 
main part of the vertical variation of ayy is not determined 
by magnetic and buoyancy effects, but rather by the kine- 
matic contribution proportional to —UxUz,y. 

The maintenance of the magnetic field over resistive 
time scales can only be explained if there is a mean electro- 
motive in the y direction that can balance the resistive losses 
of Bx. Those losses are proportional to rjdBx/dz, which 
gives the main contribution to Jy, that is proportional to 
£y{z). This is clearly suspicious and suggests that the fluc- 
tuations in the y direction are insufficient to explain the ob- 
served mean magnetic field. This result is also in agreement 
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Fig. 15 All components of the a and 77 tensors near the end of 
the run for Run A. Note that ayy and r^yx are at the noise level and 
have been multiplied by large factors (xlO^^) to see at least the 
noise. 

with the fact that ayy and r]xy, are found to be essentially 
around zero (Fig. 15). Thus, with the present knowledge we 
cannot propose a physical mechanism for the maintenance 
of the observed mean field. 

In Fig. 16 (upper row) we plot £y{z) and confirm that it 
is large enough to balance rjdBx/dz, which is overplotted 
by a dashed line. The two curves are not in perfect agree- 
ment, but this could partly be explained by the fact that 
the mean field is time-dependent and shows persistent low- 
amplitude oscillations. However, there is a caveat in that the 
definition of the fluctuating quantities, u and b that enter 
£ = u X b, have been computed using xy averages. If we 
use only y averages, so that the mean field depends on both 
X and y, we find values that are essentiaUy compatible with 
zero (see the lower row of Fig. 16). Also the strength of 
fluctuations is rather low in that case. 

7 Conclusions 

In an attempt to study the possibihty of dynamo action 
mediated by magnetic buoyancy and driven by work done 
against an externally imposed linear shear we came across a 
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number of unusual phenomena. Most surprising is probably 
the build-up and apparent maintenance of a strong magnetic 
field. The ratio of the Alfven speed to the sound speed is 
about 0.1, while the ratio of the Alfven speed to the rms 
velocity of the motions is about 10. The magnetic field is 
nearly axisymmetric, i.e. independent of y. In view of Cowl- 
ing's theorem, this clearly raises doubts about the validity 
of this result. The possibility of numerical problems is sup- 
ported by the fact that at least for Run B there is an unex- 
plained imbalance between magnetic energy input and out- 
put. On the other hand, this rather obvious imbalance ap- 
plies only to Run B, which is nearly isentropic. For Run A 
the overall energy balance is better Also, there are no obvi- 
ous deficiencies that one would normally expect in the case 
of insufficient numerical resolution such as oscillations on 
the grid scale. 

An important clue might be that the magnitude of the 
mean magnetic field in the cross-stream direction can only 
be explained if the m x 6 correlation of the fluctuations are 
evaluated as the departures from a horizontal xy average. 
If one takes just toroidal or y averages, the resulting corre- 
lation drops by six orders of magnitude into the numerical 
noise. In this context it should be noted that our calculations 
are normally performed in single precision, which has been 



sufficient for most previous applications. The present simu- 
lations may present an exception to this. 

Yet another possible clue comes from the fact that our 
models exhibit sustained oscillations. They are most promi- 
nent in the isothermal case (Run A) and less prominent in 
the nearly isentropic case (Run B). This may suggest that 
buoyancy oscillations could contribute to driving the mean 
electromotive force responsible for sustaining the mean 
field Bx in the cross-streeim direction. (The mean field in 
the streamwise direction. By, is readily explained by the 
shear.) On the other hand, an analysis of the flow in terms 
of the test-field method reveals that both ayy and rjyx are 
essentially zero in the final state, making it impossible to 
explain the magnetic field as a result of an a effect or a 
shear-current effect. 

In any case, it is clear that shear-driven flows in strati- 
fied systems can exhibit rich behavior. Magnetic buoyancy 
effects are clearly seen at early times shortly after injecting 
the initial cross-stream magnetic field. Within the present 
setup we have not seen the vigorous dynamo action with vis- 
ibly nonaxisymmetric fields reported by Cline et al. (2003). 
This could simply be related to differences in the parame- 
ters or to differences in the flow geometry. We recall that in 
Cline et al. (2003) there was a non-shearing upper part that 
is absent in the present work. In addition, they used a peri- 
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odic shear profile rather than a linear one. In any case, it will 
be worthwhile performing new simulations using perhaps 
also other methods and certainly larger resolution. It may 
be worthwhile to continue these studies using an isothermal 
setup rather than a nearly isentropic one, provided the os- 
cillations seen in the present work are indeed an important 
element of the overall dynamics. 
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